# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 13_claims_analyses.R
# Purpose: Analyses of using force in claims

# Note: must have read in data before, using 11_claims_readdata.R

# File structure:
# 1. Regression estimates, including some posterior predictive checks
# 2. Bayesian Model Averaging
# 3. IGO characteristics
# 4. Regression estimates for claims, with multiple claims per dyad-year reduced to most violent claim

# setwd("...")

library("rstan")
library("rstanarm")
library("texreg")
library("ggplot2")
library("xtable")
library("gridExtra")
library("loo")
library("lme4")
library("dplyr")
library("reshape2")
library("BMA")

# Source in functions
source("Functions/MCMClogit.fd.mat.R")
source("Functions/theme_jk.R")
source("Functions/plotBMA.R")
source("Functions/MCMC_roc_prc.r")
source("Functions/geom_flat_violin.R")

#############################################
# 1. Regression estimates, including some posterior predictive checks
#############################################

# Unit of analysis: claims

# Model 1 (for Figure 1)
m1_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat)
m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
               data = (m1_mf),
               # family = binomial(link = "logit"),
               # prior = student_t(df = 7, location = 0, scale = 2.5), 
               # prior_intercept = normal(0, 10),
               family = binomial(link = "logit"),
               prior = normal(location = 0, scale = 10), 
               prior_intercept = normal(0, 10),
               seed = 20170207,
               cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/claims_m1_stanglm.RDS")

# Table
m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:9, 1, 11), ])
m1_tab$varname <- c("IGOs with high leverage", 
    "Intangible salience of claim", 
    "Tangible salience of claim", 
    "Territorial claim", 
    "Joint democracy", 
    "Strategic rivalry",
    "UNGA ideal point difference",
    "Allies",
    "Intercept",
    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/claims_m1_table.tex", include.rownames = FALSE)

# Predicted probabilities for discussion in text
source("Functions/MCMC_simcase_probs.R")
source("Functions/MCMC_observed_probs.R")
m1_mm <- model.matrix(m1_eq, data = dat)
temp_sims <- as.matrix(m1_stanglm)
m1_simcase_probs <- MCMC_simcase_probs(model_matrix = m1_mm, 
                                       mcmc_out = temp_sims, 
                                       x_col = 2, 
                                       x_range_vec = c(0:9), 
                                       link = "logit", 
                                       lower = 0.05, upper = 0.95)
m1_observed_probs <- MCMC_observed_probs(model_matrix = m1_mm, 
                                       mcmc_out = temp_sims, 
                                       x_col = 2, 
                                       x_range_vec = c(0:9), 
                                       link = "logit", 
                                       lower = 0.05, upper = 0.95)

# First differences (Figure 1)
m1_mm <- model.matrix(m1_eq, data = dat)

temp_sims <- as.matrix(m1_stanglm)
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                           credint = c(0.025, 0.975), 
                           percentiles = c(0.1, 0.9),
                           full_sims = TRUE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]   [,2]
# [1,]     NA     NA
# [2,] 0.0000 4.0000
# [3,] 1.0000 2.3000
# [4,] 1.0000 5.0000
# [5,] 0.0000 1.0000
# [6,] 0.0000 1.0000
# [7,] 0.0000 1.0000
# [8,] 0.1496 2.6827
# [9,] 0.0000 1.0000

fd.dat <- data.frame(m1.fd)
fd.dat$variable_label <- c("IGOs with high leverage\n(from 0 to 4)", 
                           "Intangible salience of claim\n(from 1 to 2.3)", 
                           "Tangible salience of claim\n(from 1 to 5)", 
                           "Territorial claim\n(no to yes)", 
                           "Joint democracy\n(no to yes)", 
                           "Strategic rivalry\n(no to yes)",
                           "UNGA ideal point difference\n(0.1 to 2.7)",
                           "Allies \n(no to yes)")


fd_long <- melt(m1.fd)

fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", 1,
                         ifelse(fd_long$Var2 == "salint", 2, 
                                ifelse(fd_long$Var2 == "saltan", 3,
                                       ifelse(fd_long$Var2 == "terriss", 4,
                                              ifelse(fd_long$Var2 == "jointpol7", 5,
                                                     ifelse(fd_long$Var2 == "rivalry_th", 6, 
                                                            ifelse(fd_long$Var2 == "absidealdiff", 7,
                                                                   ifelse(fd_long$Var2 == "atopally", 8, NA))))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", "IGOs with high leverage\n(from 0 to 4)",
                                 ifelse(fd_long$Var2 == "salint", "Intangible salience of claim\n(from 1 to 2.3)", 
                                        ifelse(fd_long$Var2 == "saltan", "Tangible salience of claim\n(from 1 to 5)",
                                               ifelse(fd_long$Var2 == "terriss", "Territorial claim\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)", 
                                                                    ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 2.7)",
                                                                           ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)", NA))))))))

# ROPE: use SE of the ratio of 1s to 0s
r_n <- length(m1_mf$useforce)
r_p <- sum(m1_mf$useforce) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of using force\nduring a claim as each predictor changes as indicated") + 
  xlab("") + 
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Use of force less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Use of force more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.8, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-3, 3]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = ROPE, x = 7.25, xend = 7.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 8.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 8.25, xend = 8.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 3.5, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 3.45, xend = 3.15, size = 0.25) +
  annotate(geom = "text", y = -0.5, x = 7.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.275, x = 7.65, xend = 7.9, size = 0.25)

# % < 0
m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.1, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p, file = "Output_Tables-and-Figures/claims_m1_fd.pdf", width = 6, height = 8)

# Remove transparency
fd.p2 <- fd.p + annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray90", color = NA) + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.1, size = 3)  + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)
ggsave(fd.p2, file = "Output_Tables-and-Figures/claims_m1_fd.eps", width = 6, height = 8)


# LOO
m1_noigo <- stan_glm(formula = useforce ~ salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally,
                     data = (m1_mf),
                     family = binomial(link = "logit"),
                     prior = normal(0, 10), 
                     prior_intercept = normal(0, 10),
                     cores = 4,
                     seed = 20170207)

saveRDS(m1_noigo, file = "Output_MCMC/claims_m1_noigo_stanglm.RDS")

loo_noigo <- loo(m1_noigo, cores = 2)
loo_m1 <- loo(m1_stanglm, cores = 2)
compare(loo_noigo, loo_m1)
# elpd_diff        se 
# 1.0       2.1 

# Influence of individual observations
p <- ggplot(data = data.frame(x = 1:168, y = loo_m1$pareto_k),
            aes(x = x, y = y)) + geom_point(shape = 1) + xlab("Observation") + ylab("Shape parameter k") + geom_hline(yintercept = 0.5, linetype = "dashed") + theme_jk()
ggsave(p, file = "Output_Tables-and-Figures/claims_m1_loo_k.pdf", width = 5, height = 4)

# Other model fit diagnostics via Stan

# Data 
roc_mf <-  model.frame(useforce ~ igo_lev3_count_use_bc + igo_ingr23_use_bc + igo_mp3_use_bc + igo_pb_use + igo_allothers_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + name + claim, data = dat)
N <- nrow(roc_mf)
useforce <- roc_mf$useforce
igo_lev3_count_use_bc <- roc_mf$igo_lev3_count_use_bc
igo_ingr23_use_bc <- roc_mf$igo_ingr23_use_bc
igo_mp3_use_bc <- roc_mf$igo_mp3_use_bc
igo_pb_use <- roc_mf$igo_pb_use
igo_allothers_bc <- roc_mf$igo_allothers_bc
salint <- roc_mf$salint
saltan <- roc_mf$saltan
terriss <- roc_mf$terriss
jointpol7 <- roc_mf$jointpol7
rivalry_th <- roc_mf$rivalry_th
absidealdiff <- roc_mf$absidealdiff
atopally <- roc_mf$atopally

m1_roc_data_list <- c("N", "useforce", "igo_lev3_count_use_bc", "igo_ingr23_use_bc", "igo_mp3_use_bc", "igo_pb_use", "igo_allothers_bc", "salint", "saltan", "terriss", "jointpol7", "rivalry_th", "absidealdiff", "atopally")

# Fit m1 directly in Stan
m1_roc_rstan <- stan(file = "claims_m1_roc.stan", 
                   data = m1_roc_data_list,
                   iter = 1000, chains = 4,
                   cores = 4,
                   seed = 20170207)

print(m1_roc_rstan)
saveRDS(m1_roc_rstan, file = "Output_MCMC/claims_m1_rstan.RDS")

# Fit igo_ingr23_use_bc directly in Stan
m_ingr23_roc_rstan <- stan(file = "claims_m-ingr23_roc.stan", 
                     data = m1_roc_data_list,
                     iter = 1000, chains = 4,
                     seed = 20170207)
saveRDS(m_ingr23_roc_rstan, file = "Output_MCMC/claims_m1_ingr23_rstan.RDS")

# Fit igo_mp3_use_bc directly in Stan
m_mp3_roc_rstan <- stan(file = "claims_m-mp3_roc.stan", 
                     data = m1_roc_data_list,
                     iter = 1000, chains = 4,
                     seed = 20170207)
saveRDS(m_mp3_roc_rstan, file = "Output_MCMC/claims_m1_mp3_rstan.RDS")

# Fit igo_pb_use directly in Stan
m_pb_roc_rstan <- stan(file = "claims_m-pb_roc.stan", 
                     data = m1_roc_data_list,
                     iter = 1000, chains = 4,
                     seed = 20170207)
saveRDS(m_pb_roc_rstan, file = "Output_MCMC/claims_m1_pb_rstan.RDS")

# Fit all other IGOs directly in Stan
m_allothers_roc_rstan <- stan(file = "claims_m-allothers_roc.stan", 
                              data = m1_roc_data_list,
                              iter = 1000, chains = 4,
                              seed = 20170207)
saveRDS(m_allothers_roc_rstan, file = "Output_MCMC/claims_m1_allothers_rstan.RDS")

# Fit empty model (no IGOs) directly in Stan
m_noigos_roc_rstan <- stan(file = "claims_m-noigos_roc.stan", 
                           data = m1_roc_data_list,
                           iter = 1000, chains = 4,
                           seed = 20170207)
saveRDS(m_noigos_roc_rstan, file = "Output_MCMC/claims_m1_noigos_rstan.RDS")

# ROC and precision-recall curves

m1_rp <- MCMC_roc_prc(stan_object = m1_roc_rstan, 
                      model_frame = m1_mf, 
                      model_name = "IGOs with high leverage")

m_ingr23_rp <- MCMC_roc_prc(stan_object = m_ingr23_roc_rstan, 
                               model_frame = m1_mf, 
                               model_name = "Structured IGOs")

m_mp3_rp <- MCMC_roc_prc(stan_object = m_mp3_roc_rstan, 
                               model_frame = m1_mf, 
                               model_name = "Highly structured Security IGOs")

m_pb_rp <- MCMC_roc_prc(stan_object = m_pb_roc_rstan, 
                               model_frame = m1_mf, 
                               model_name = "Peace-brokering IGOs")

m_allothers_rp <- MCMC_roc_prc(stan_object = m_allothers_roc_rstan, 
                    model_frame = m1_mf, 
                    model_name = "All other IGOs")

m_noigos_rp <- MCMC_roc_prc(stan_object = m_noigos_roc_rstan, 
                      model_frame = m1_mf, 
                      model_name = "No IGOs")

PRC_dat <- rbind(m1_rp$prc_dat, m_ingr23_rp$prc_dat, m_mp3_rp$prc_dat, m_pb_rp$prc_dat, m_allothers_rp$prc_dat, m_noigos_rp$prc_dat)
PRC_labeldat <- rbind(m1_rp$area_under_prc, m_ingr23_rp$area_under_prc, m_mp3_rp$area_under_prc, m_pb_rp$area_under_prc, m_allothers_rp$area_under_prc, m_noigos_rp$area_under_prc)

PRC.p <- ggplot(data = PRC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = PRC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 0, y = 0.2, hjust = 0) + labs(title = "Precision-recall curves", hjust = 0.5) + xlab("Recall") + ylab("Precision") + theme_jk()

ROC_dat <- rbind(m1_rp$roc_dat, m_ingr23_rp$roc_dat, m_mp3_rp$roc_dat, m_pb_rp$roc_dat, m_allothers_rp$roc_dat, m_noigos_rp$roc_dat)
ROC_labeldat <- rbind(m1_rp$area_under_roc, m_ingr23_rp$area_under_roc, m_mp3_rp$area_under_roc, m_pb_rp$area_under_roc, m_allothers_rp$area_under_roc, m_noigos_rp$area_under_roc)

ROC.p <- ggplot(data = ROC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = ROC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 1, y = 0.2, hjust = 1) + labs(title = "ROC curves", hjust = 0.5) + xlab("1 - Specificity") + ylab("Sensitivity") + theme_jk()

pdf(file = "Output_Tables-and-Figures/claims_pr-roc.pdf", width = 6, height = 12)
grid.arrange(PRC.p, ROC.p, ncol = 2)
dev.off()

# Posterior predictive checks

m1_stan_object <- m1_roc_rstan
m1_pred_sim <- as.data.frame(m1_stan_object, pars = c("pred"))
m1_pred_sim <- reshape2::melt(m1_pred_sim)
m1_pred_sim$obs <- as.numeric(gsub("[^0-9]","",as.character(m1_pred_sim$variable)))

m1_obs_igo <- data.frame(y_obs = m1_mf[, "useforce"], 
                      x_obs = m1_mf[, "igo_lev3_count_use_bc"],
                      obs = 1:(nrow(m1_mf)))

m1_median_obs <- summarize(group_by(m1_obs_igo, x_obs),
                        median_y = mean(y_obs))

m1_pp_df <- merge(x = m1_pred_sim[, c("value", "obs")], y = m1_obs_igo, by = "obs")

ppc.p <- ggplot(data = m1_pp_df, aes(x = x_obs, y = value, group = x_obs)) + geom_boxplot(outlier.size = NA, size = 0.25) + scale_x_continuous(breaks = unique(m1_obs_igo$x_obs)) + geom_point(data = m1_median_obs, aes(x = x_obs, y = median_y), color = "red") + theme_jk() + xlab("Joint memberships in IGOs with high leverage") + ylab("Estimated probability of using force") + annotate("point", x = 0.2, y = 0.9, color = "red") + annotate("text", x = 0.4, y = 0.9, hjust = 0, label = "Percent of observed claims\nexperiencing the use of force", size = 3) + annotate("boxplot", x = 5, lower = 0.85, upper = 0.95, middle = 0.9, ymin = 0.8, ymax = 1, size = 0.25) + annotate("text", x = 5.5, y = 0.9, hjust = 0, label = "Posterior predictive\ndistribution", size = 3)

ggsave(ppc.p, file = "Output_Tables-and-Figures/claims_ppc.pdf", width = 6, height = 4)

# Models with other controls

# igo_ingr23_use_bc
m2_eq <- useforce ~ igo_lev3_count_use_bc + igo_ingr23_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m2_mf <- model.frame(m2_eq, data = dat)
m2 <- glm(m2_eq, data = m2_mf, family = binomial(link = "logit"))

m2_stanglm <- stan_glm(formula = m2_eq,
                       data = (m2_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m2_stanglm, file = "Output_MCMC/claims_m2_stanglm.RDS")

# Table
m2_tab <- m2_stanglm$stan_summary[, c(1, 3)]
m2_tab <- data.frame(m2_tab[c(2:10, 1, 12), ])
m2_tab$varname <- c("IGOs with high leverage",
                    "Structured IGOs", 
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m2_tab) <- c("Mean", "Standard deviation", "Variable")
m2_tab <- m2_tab[, c(3, 1, 2)]

m2_xtab <- xtable(m2_tab, digits = 2)
print(m2_xtab, file = "Output_MCMC/claims_m2_table.tex", include.rownames = FALSE)

# with HLIGOs explicitly removed
m2b_eq <- useforce ~ igo_lev3_count_use_bc + I(igo_ingr23_use_bc - igo_lev3_count_use_bc) + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m2b_mf <- model.frame(m2b_eq, data = dat)
m2b <- glm(useforce ~ ., data = m2b_mf, family = binomial(link = "logit"))

texreg(list(m2b),
       stars = 0.1,
       file = "Output_Tables-and-Figures/claims_m2b_mle_table.tex",
       custom.note = "* p < 0.05, one-tailed tests.",
       digits = 2,
       caption = "Determinants of using force in claims: comparing IGOs with high leverage to Structured IGOs. Logistic regression estimates (fit with maximum likelihood)",
       caption.above = TRUE,
       custom.coef.names = c("Intercept",
                             "IGOs with high leverage",
                             "Structured IGOs (without IGOs w/ high leverage)",
                             "Intangible salience of claim", 
                             "Tangible salience of claim", 
                             "Territorial claim", 
                             "Joint democracy", 
                             "Strategic rivalry",
                             "UNGA ideal point difference",
                             "Allies"))

# igo_mp3_use_bc
m3_eq <- useforce ~ igo_lev3_count_use_bc + igo_mp3_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m3_mf <- model.frame(m3_eq, data = dat)
m3 <- glm(m3_eq, data = m3_mf, family = binomial(link = "logit"))

m3_stanglm <- stan_glm(formula = m3_eq,
                       data = (m3_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m3_stanglm, file = "Output_MCMC/claims_m3_stanglm.RDS")

# Table
m3_tab <- m3_stanglm$stan_summary[, c(1, 3)]
m3_tab <- data.frame(m3_tab[c(2:10, 1, 12), ])
m3_tab$varname <- c("IGOs with high leverage",
                    "Highly structured Security IGOs", 
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m3_tab) <- c("Mean", "Standard deviation", "Variable")
m3_tab <- m3_tab[, c(3, 1, 2)]

m3_xtab <- xtable(m3_tab, digits = 2)
print(m3_xtab, file = "Output_MCMC/claims_m3_table.tex", include.rownames = FALSE)

# igo_pb_use
m4_eq <- useforce ~ igo_lev3_count_use_bc + igo_pb_use + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m4_mf <- model.frame(m4_eq, data = dat)
m4 <- glm(m4_eq, data = m4_mf, family = binomial(link = "logit"))

m4_stanglm <- stan_glm(formula = m4_eq,
                       data = (m4_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m4_stanglm, file = "Output_MCMC/claims_m4_stanglm.RDS")

# Table
m4_tab <- m4_stanglm$stan_summary[, c(1, 3)]
m4_tab <- data.frame(m4_tab[c(2:10, 1, 12), ])
m4_tab$varname <- c("IGOs with high leverage",
                    "Peace-brokering IGOs", 
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m4_tab) <- c("Mean", "Standard deviation", "Variable")
m4_tab <- m4_tab[, c(3, 1, 2)]

m4_xtab <- xtable(m4_tab, digits = 2)
print(m4_xtab, file = "Output_Tables-and-Figures/claims_m4_table.tex", include.rownames = FALSE)

# igo_allothers_bc
m5_eq <- useforce ~ igo_lev3_count_use_bc + igo_allothers_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m5_mf <- model.frame(m5_eq, data = dat_spec)
m5 <- glm(m5_eq, data = m5_mf, family = binomial(link = "logit"))

m5_stanglm <- stan_glm(formula = m5_eq,
                       data = (m5_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m5_stanglm, file = "Output_MCMC/claims_m5_stanglm.RDS")

# Table
m5_tab <- m5_stanglm$stan_summary[, c(1, 3)]
m5_tab <- data.frame(m5_tab[c(2:10, 1, 12), ])
m5_tab$varname <- c("IGOs with high leverage",
                    "All other IGOs", 
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m5_tab) <- c("Mean", "Standard deviation", "Variable")
m5_tab <- m5_tab[, c(3, 1, 2)]

m5_xtab <- xtable(m5_tab, digits = 2)
print(m5_xtab, file = "Output_Tables-and-Figures/claims_m5_table.tex", include.rownames = FALSE)

# cincdif
m6_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + cincdif
m6_mf <- model.frame(m6_eq, data = dat)
m6 <- glm(m6_eq, data = m6_mf, family = binomial(link = "logit"))

m6_stanglm <- stan_glm(formula = m6_eq,
                       data = (m6_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m6_stanglm, file = "Output_MCMC/claims_m6_stanglm.RDS")

# Table
m6_tab <- m6_stanglm$stan_summary[, c(1, 3)]
m6_tab <- data.frame(m6_tab[c(2:10, 1, 12), ])
m6_tab$varname <- c("IGOs with high leverage",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Power differential",
                    "Intercept",
                    "Log-posterior density")
names(m6_tab) <- c("Mean", "Standard deviation", "Variable")
m6_tab <- m6_tab[, c(3, 1, 2)]

m6_xtab <- xtable(m6_tab, digits = 2)
print(m6_xtab, file = "Output_Tables-and-Figures/claims_m6_table.tex", include.rownames = FALSE)

# tradedep_low_ln
m7_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + tradedep_low_ln
m7_mf <- model.frame(m7_eq, data = dat)
m7 <- glm(m7_eq, data = m7_mf, family = binomial(link = "logit"))

m7_stanglm <- stan_glm(formula = m7_eq,
                       data = (m7_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m7_stanglm, file = "Output_MCMC/claims_m7_stanglm.RDS")

# Table
m7_tab <- m7_stanglm$stan_summary[, c(1, 3)]
m7_tab <- data.frame(m7_tab[c(2:10, 1, 12), ])
m7_tab$varname <- c("IGOs with high leverage",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Trade dependence (lower)",
                    "Intercept",
                    "Log-posterior density")
names(m7_tab) <- c("Mean", "Standard deviation", "Variable")
m7_tab <- m7_tab[, c(3, 1, 2)]

m7_xtab <- xtable(m7_tab, digits = 2)
print(m7_xtab, file = "Output_Tables-and-Figures/claims_m7_table.tex", include.rownames = FALSE)

# gdppc_low_ln
m8_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + gdppc_low_ln
m8_mf <- model.frame(m8_eq, data = dat)
m8 <- glm(m8_eq, data = m8_mf, family = binomial(link = "logit"))

m8_stanglm <- stan_glm(formula = m8_eq,
                       data = (m8_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m8_stanglm, file = "Output_MCMC/claims_m8_stanglm.RDS")

# Table
m8_tab <- m8_stanglm$stan_summary[, c(1, 3)]
m8_tab <- data.frame(m8_tab[c(2:10, 1, 12), ])
m8_tab$varname <- c("IGOs with high leverage",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "GDP per capita (lower)",
                    "Intercept",
                    "Log-posterior density")
names(m8_tab) <- c("Mean", "Standard deviation", "Variable")
m8_tab <- m8_tab[, c(3, 1, 2)]

m8_xtab <- xtable(m8_tab, digits = 2)
print(m8_xtab, file = "Output_Tables-and-Figures/claims_m8_table.tex", include.rownames = FALSE)

# regionFactor
m9_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + regionFactor
m9_mf <- model.frame(m9_eq, data = dat)
m9 <- glm(m9_eq, data = m9_mf, family = binomial(link = "logit"))

m9_stanglm <- stan_glm(formula = m9_eq,
                       data = (m9_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m9_stanglm, file = "Output_MCMC/claims_m9_stanglm.RDS")

# Table
m9_tab <- m9_stanglm$stan_summary[, c(1, 3)]
m9_tab <- data.frame(m9_tab[c(2:11, 1, 13), ])
m9_tab$varname <- c("IGOs with high leverage",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Europe (vs. Western Hemisphere)",
                    "Middle East (vs. Western Hemisphere)",
                    "Intercept",
                    "Log-posterior density")
names(m9_tab) <- c("Mean", "Standard deviation", "Variable")
m9_tab <- m9_tab[, c(3, 1, 2)]

m9_xtab <- xtable(m9_tab, digits = 2)
print(m9_xtab, file = "Output_Tables-and-Figures/claims_m9_table.tex", include.rownames = FALSE)

# coldWar
m10_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + coldWar
m10_mf <- model.frame(m10_eq, data = dat)
m10 <- glm(m10_eq, data = m10_mf, family = binomial(link = "logit"))

m10_stanglm <- stan_glm(formula = m10_eq,
                        data = (m10_mf),
                        family = binomial(link = "logit"),
                        prior = normal(location = 0, scale = 10), 
                        prior_intercept = normal(0, 10),
                        seed = 20170207,
                        cores = 4)

saveRDS(m10_stanglm, file = "Output_MCMC/claims_m10_stanglm.RDS")

# Table
m10_tab <- m10_stanglm$stan_summary[, c(1, 3)]
m10_tab <- data.frame(m10_tab[c(2:10, 1, 12), ])
m10_tab$varname <- c("IGOs with high leverage",
                     "Intangible salience of claim", 
                     "Tangible salience of claim", 
                     "Territorial claim", 
                     "Joint democracy", 
                     "Strategic rivalry",
                     "UNGA ideal point difference",
                     "Allies",
                     "Cold War",
                     "Intercept",
                     "Log-posterior density")
names(m10_tab) <- c("Mean", "Standard deviation", "Variable")
m10_tab <- m10_tab[, c(3, 1, 2)]

m10_xtab <- xtable(m10_tab, digits = 2)
print(m10_xtab, file = "Output_Tables-and-Figures/claims_m10_table.tex", include.rownames = FALSE)

# Old salience index
m14_eq <- useforce ~ igo_lev3_count_use_bc + salold + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m14_mf <- model.frame(m14_eq, data = dat)
m14 <- glm(m14_eq, data = m14_mf, family = binomial(link = "logit"))

m14_stanglm <- stan_glm(formula = m14_eq,
                        data = (m14_mf),
                        family = binomial(link = "logit"),
                        prior = normal(location = 0, scale = 10), 
                        prior_intercept = normal(0, 10),
                        seed = 20170207,
                        cores = 4)

saveRDS(m14_stanglm, file = "Output_MCMC/claims_m14_stanglm.RDS")

# Table
m14_tab <- m14_stanglm$stan_summary[, c(1, 3)]
m14_tab <- data.frame(m14_tab[c(2:8, 1, 10), ])
m14_tab$varname <- c("IGOs with high leverage",
                     "Salience of claim (old index)", 
                     "Territorial claim", 
                     "Joint democracy", 
                     "Strategic rivalry",
                     "UNGA ideal point difference",
                     "Allies",
                     "Intercept",
                     "Log-posterior density")
names(m14_tab) <- c("Mean", "Standard deviation", "Variable")
m14_tab <- m14_tab[, c(3, 1, 2)]

m14_xtab <- xtable(m14_tab, digits = 2)
print(m14_xtab, file = "Output_Tables-and-Figures/claims_m14_table.tex", include.rownames = FALSE)

# Random effects for time
m15_eq <- useforce ~ (igo_lev3_count_use_bc | decade) + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m15_mf <- model.frame(useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + decade, data = dat)
m15 <- glmer(m15_eq, data = dat, family = binomial(link = "logit"))

m15_stanglm <- stan_glmer(formula = m15_eq,
                        data = (dat),
                        family = binomial(link = "logit"),
                        prior = normal(location = 0, scale = 10), 
                        prior_intercept = normal(0, 10),
                        seed = 20170207,
                        cores = 4)

saveRDS(m15_stanglm, file = "Output_MCMC/claims_m15_stanglm.RDS")

m15_mcmc <- as.data.frame(m15_stanglm$stanfit)
m15_re_mcmc <- m15_mcmc[ , grepl( "igo_lev3_count_use_bc decade:" , names(m15_mcmc) ) ]
m15_re_mcmc$`b[igo_lev3_count_use_bc decade:_NEW_decade]` <- NULL
table(m15_mf$decade)
# 0  1  2  3  4  5 
# 6 31 27 35 31 38 
names(m15_re_mcmc) <- c("1947-1949 (6 claims)", "1950-1959 (31 claims)", "1960-1969 (27 claims)", "1970-1979 (35 claims)", "1980-1989 (31 claims)", "1990-2001 (38 claims)")

m15_re_mcmc_probs <- data.frame(variable = names(m15_re_mcmc), 
                                perc_below_0 = apply(m15_re_mcmc, MARGIN = 2, FUN = function(x) sum(ifelse(x < 0, 1, 0)) / length(x)),
                                perc_larger = apply(m15_re_mcmc, MARGIN = 2, FUN = function(x) ifelse(sum(ifelse(x < 0, 1, 0)) / length(x) > 0.5, sum(ifelse(x < 0, 1, 0)) / length(x), sum(ifelse(x > 0, 1, 0)) / length(x))),
                                perc_below_0_label = paste(round(apply(m15_re_mcmc, MARGIN = 2, FUN = function(x) sum(ifelse(x < 0, 1, 0)) / length(x)) * 100, digits = 1), "% < 0", sep = ""), 
                                xpos = apply(m15_re_mcmc, MARGIN = 2, FUN = function(x) quantile(x, probs = 0.01)))

m15_re_mcmc_long <- reshape2::melt(m15_re_mcmc)
p <- ggplot(data = m15_re_mcmc_long, aes(x = value)) + geom_density(color = NA, fill = "lightgray") + geom_text(data = m15_re_mcmc_probs, aes(x = xpos, y = 1, label = perc_below_0_label)) + facet_wrap(~ variable, scales = "free", ncol = 1) + geom_vline(xintercept = 0) + theme_jk() + xlab("Estimate for decade-specific logit coefficients on IGOs with high leverage") + ylab("") + scale_y_continuous(breaks = NULL)

ggsave(p, file = "Output_Tables-and-Figures/claims_hligo_re.pdf", width = 6, height = 8)

# Model comparison: do the REs improve model fit? 
loo_m15 <- loo(m15_stanglm, cores = 2)
compare(loo_m1, loo_m15)
# elpd_diff   se 
# 1.9         3.4 

# Interact HLIGOs and commitment problem (power differential)
m16_eq <- useforce ~ igo_lev3_count_use_bc * cincdif + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m16_mf <- model.frame(m16_eq, data = dat)
m16 <- glm(m16_eq, data = m16_mf, family = binomial(link = "logit"))

m16_mf$cincdif_q3 <- cut(m16_mf$cincdif, breaks = quantile(m16_mf$cincdif, probs = c(0, 0.33, 0.66, 1)), include.lowest = TRUE)
levels(m16_mf$cincdif_q3) <- c("[0, 0.003]", "(0.003, 0.03]", "(0.03, 0.31]")

m16b_eq <- useforce ~ igo_lev3_count_use_bc * cincdif_q3 + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m16b <- glm(m16b_eq, data = m16_mf, family = binomial(link = "logit"))


m16_stanglm <- stan_glm(formula = m16b_eq,
                       data = (m16_mf),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m16_stanglm, file = "Output_MCMC/claims_m16_stanglm.RDS")

m16_mcmc <- as.data.frame(m16_stanglm$stanfit)
igo_lowcincdif <- data.frame(x = m16_mcmc[, "igo_lev3_count_use_bc"], var = "No or small CINC difference (0 to 0.003)")
igo_midcincdif <- data.frame(x = m16_mcmc[, "igo_lev3_count_use_bc"] + m16_mcmc[, "igo_lev3_count_use_bc:cincdif_q3(0.003, 0.03]"], var = "Medium CINC difference (0.003 to 0.03)")
igo_hicincdif <- data.frame(x = m16_mcmc[, "igo_lev3_count_use_bc"] + m16_mcmc[, "igo_lev3_count_use_bc:cincdif_q3(0.03, 0.31]"], var = "High CINC difference (0.03 to 0.31)")
m16_interact_mcmc <- rbind(igo_lowcincdif, igo_midcincdif, igo_hicincdif)

perc_below_0 <- c(sum(ifelse(igo_lowcincdif$x < 0, 1, 0)) / length(igo_lowcincdif$x), sum(ifelse(igo_midcincdif$x < 0, 1, 0)) / length(igo_midcincdif$x), sum(ifelse(igo_hicincdif$x < 0, 1, 0)) / length(igo_hicincdif$x))

m16_interact_mcmc_probs <- data.frame(var = c("No or small CINC difference (0 to 0.003)", "Medium CINC difference (0.003 to 0.03)", "High CINC difference (0.03 to 0.31)"), 
                                perc_below_0 = c(sum(ifelse(igo_lowcincdif$x < 0, 1, 0)) / length(igo_lowcincdif$x), sum(ifelse(igo_midcincdif$x < 0, 1, 0)) / length(igo_midcincdif$x), sum(ifelse(igo_hicincdif$x < 0, 1, 0)) / length(igo_hicincdif$x)),
                                perc_below_0_label = paste(round(perc_below_0 * 100, digits = 1), "% < 0", sep = ""),
                                xpos = c(-1, -1, -1.5))

p <- ggplot(data = m16_interact_mcmc, aes(x = x)) + geom_density(color = NA, fill = "lightgray") + geom_text(data = m16_interact_mcmc_probs, aes(x = xpos, y = 0.6, label = perc_below_0_label)) + facet_wrap(~ var, scales = "free", ncol = 1) + geom_vline(xintercept = 0) + theme_jk() + xlab("Estimate for conditional effect (dy/dx) of IGOs with high leverage") + ylab("") + scale_y_continuous(breaks = NULL)

ggsave(p, file = "Output_Tables-and-Figures/claims_igo-cinc_interact.pdf", width = 5, height = 8)

# Table of MLE results
texreg(list(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m14),
          stars = 0.1,
        file = "Output_Tables-and-Figures/claims_mle-results.tex",
        custom.note = "* p < 0.05, one-tailed tests.",
        digits = 2,
        caption = "Determinants of using force in claims: logistic regression estimates (fit with maximum likelihood)",
        caption.above = TRUE,
        custom.coef.names = c("Intercept",
                              "IGOs with high leverage",
                              "Intangible salience of claim", 
                              "Tangible salience of claim", 
                              "Territorial claim", 
                              "Joint democracy", 
                              "Strategic rivalry",
                              "UNGA ideal point difference",
                              "Allies",
                              "Structured IGOs",
                              "Highly structured Security IGOs", 
                              "Peace-brokering IGOs",
                              "All other IGOs",
                              "Power differential",
                              "Trade dependence (lower)",
                              "GDP per capita (lower)",
                              "Europe (vs. Western Hemisphere)",
                              "Middle East (vs. Western Hemisphere)",
                              "Cold war",
                              "Salience of claim (aggregate measure)"))

#############################################
# 2. Bayesian Model Averaging
#############################################

bma_dat <- model.frame(useforce ~ igo_lev3_count_use_bc + igo_ingr23_use_bc + igo_mp3_use_bc + igo_pb_use + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + cincdif, data = dat)
names(bma_dat) <- c("Use of force",
                    "IGOs with high leverage",
                    "Structured IGOs",
                    "Highly structured Security IGOs", 
                    "Peace-brokering IGOs",
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Power differential")

claims_bma <- bic.glm(y = bma_dat[, 1], x = bma_dat[, c(2:13)], 
                      glm.family = binomial(link = "logit"))

pdf("Output_Tables-and-Figures/claims_bma.pdf", width = 12, height = 8)
plotJK.bic.glm(claims_bma, mfrow = c(3, 4))
dev.off()

#############################################
# Unit of analysis: claim-years: See claims_dy_analyses.R
############################################## 

#############################################
# 3. IGO characteristics
#############################################

dat$idealpointdiff_std <- scales::rescale(x = dat$absidealdiff, to = c(0, 1))
dat$cincdif_std <- scales::rescale(x = dat$cincdif, to = c(0, 1))

hligo_sum <- summarize(group_by(dat[is.na(dat$igo_lev3_count_use_bc) == FALSE, ], igo_lev3_count_use_bc), 
                       obs_count = length(atopally) / nrow(dat[is.na(dat$igo_lev3_count_use_bc) == FALSE, ]), 
                       democracies_mean = mean(jointpol7, na.rm = TRUE), 
                       rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                       idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                       allies_mean = mean(atopally, na.rm = TRUE), 
                       cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(hligo_sum) <- c("Joint IGO memberships", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
hligo_sum$`IGO type` <- 1 # "IGOs with high leverage"

# igo_ingr23_use_bc
ingr23_sum <- summarize(group_by(dat[is.na(dat$igo_ingr23_use_bc) == FALSE, ], igo_ingr23_use_bc), 
                        obs_count = length(atopally) / nrow(dat[is.na(dat$igo_lev3_count_use_bc) == FALSE, ]), 
                        democracies_mean = mean(jointpol7, na.rm = TRUE), 
                        rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                        idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                        allies_mean = mean(atopally, na.rm = TRUE), 
                        cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(ingr23_sum) <- c("Joint IGO memberships", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
ingr23_sum$`IGO type` <- 2 # "Structured IGOs"

# igo_mp3_use_bc
mp3_sum <- summarize(group_by(dat[is.na(dat$igo_mp3_use_bc) == FALSE, ], igo_mp3_use_bc), 
                     obs_count = length(atopally) / nrow(dat[is.na(dat$igo_lev3_count_use_bc) == FALSE, ]), 
                     democracies_mean = mean(jointpol7, na.rm = TRUE), 
                     rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                     idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                     allies_mean = mean(atopally, na.rm = TRUE), 
                     cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(mp3_sum) <- c("Joint IGO memberships", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
mp3_sum$`IGO type` <- 3 # "Highly structured Security IGOs"

# igo_pb_use
pb_sum <- summarize(group_by(dat[is.na(dat$igo_pb_use) == FALSE, ], igo_pb_use), 
                    obs_count = length(atopally) / nrow(dat[is.na(dat$igo_lev3_count_use_bc) == FALSE, ]), 
                    democracies_mean = mean(jointpol7, na.rm = TRUE), 
                    rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                    idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                    allies_mean = mean(atopally, na.rm = TRUE), 
                    cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(pb_sum) <- c("Joint IGO memberships", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
pb_sum$`IGO type` <- 4 # "Peace-brokering IGOs"

igo_sum <- rbind(hligo_sum, ingr23_sum, mp3_sum, pb_sum)
names(igo_sum) <- c("Joint IGO memberships", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential", "IGO type")
igo_sum$`IGO type` <- factor(igo_sum$`IGO type`, labels = c("IGOs with high leverage", "Structured IGOs", "Highly structured Security IGOs", "Peace-brokering IGOs"))
igo_sum <- igo_sum[is.na(igo_sum$`Joint IGO memberships`) == FALSE, ]

igo_sum_long <- reshape2::melt(igo_sum, id.vars = c("Joint IGO memberships", "IGO type", "Observations"))

p <- ggplot(data = igo_sum_long, aes(x = `Joint IGO memberships`, y = value)) + facet_grid(variable ~ `IGO type`, scales = "free_x") + geom_point(aes(size = Observations)) + geom_line(aes(group = 1))

# Characteristics across IGO memberships (quartiles)

dat$hligo_quartiles <- cut(dat$igo_lev3_count_use_bc, breaks = quantile(dat$igo_lev3_count_use_bc, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), include.lowest = TRUE)
dat$hligo_quartiles <- factor(dat$hligo_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

dat$ingr23_quartiles <- cut(dat$igo_ingr23_use_bc, breaks = quantile(dat$igo_ingr23_use_bc, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE), include.lowest = TRUE)
dat$ingr23_quartiles <- factor(dat$ingr23_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

# dat$mp3_quartiles <- cut(dat$igo_mp3_use_bc, breaks = quantile(dat$igo_mp3_use_bc, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE))
dat$mp3_quartiles <- ifelse(dat$igo_mp3_use_bc >= 3, 4, dat$igo_mp3_use_bc + 1)
dat$mp3_quartiles <- factor(dat$mp3_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

# dat$pb_quartiles <- cut(dat$igo_pb_use, breaks = quantile(dat$igo_pb_use, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE))
dat$pb_quartiles <- ifelse(dat$igo_pb_use <= 1, 1, 
                           ifelse(dat$igo_pb_use == 2, 2,
                                  ifelse(dat$igo_pb_use == 3, 3,
                                         ifelse(dat$igo_pb_use >= 4, 4, NA))))
dat$pb_quartiles <- factor(dat$pb_quartiles, labels = c("1st", "2nd", "3rd", "4th"))

dat$idealpointdiff_std <- scales::rescale(x = dat$idealpointdiff_std, to = c(0, 1))
dat$cincdif_std <- scales::rescale(x = dat$cincdif, to = c(0, 1))

hligo_sum <- summarize(group_by(dat[is.na(dat$hligo_quartiles) == FALSE, ], hligo_quartiles), 
                       obs_ratio = length(atopally) / nrow(dat[is.na(dat$hligo_quartiles) == FALSE, ]), 
                       obs_count = length(atopally),
                       democracies_mean = mean(jointpol7, na.rm = TRUE), 
                       rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                       idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                       allies_mean = mean(atopally, na.rm = TRUE), 
                       cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(hligo_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
hligo_sum$`IGO type` <- 1 # "IGOs with high leverage"

# igo_ingr23_use_bc
ingr23_sum <- summarize(group_by(dat[is.na(dat$ingr23_quartiles) == FALSE, ], ingr23_quartiles), 
                        obs_ratio = length(atopally) / nrow(dat[is.na(dat$ingr23_quartiles) == FALSE, ]), 
                        obs_count = length(atopally),
                        democracies_mean = mean(jointpol7, na.rm = TRUE), 
                        rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                        idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                        allies_mean = mean(atopally, na.rm = TRUE), 
                        cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(ingr23_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
ingr23_sum$`IGO type` <- 2 # "Structured IGOs"

# igo_mp3_use_bc
mp3_sum <- summarize(group_by(dat[is.na(dat$mp3_quartiles) == FALSE, ], mp3_quartiles), 
                     obs_ratio = length(atopally) / nrow(dat[is.na(dat$mp3_quartiles) == FALSE, ]), 
                     obs_count = length(atopally),
                     democracies_mean = mean(jointpol7, na.rm = TRUE), 
                     rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                     idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                     allies_mean = mean(atopally, na.rm = TRUE), 
                     cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(mp3_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
mp3_sum$`IGO type` <- 3 # "Highly structured Security IGOs"

# igo_pb_use
pb_sum <- summarize(group_by(dat[is.na(dat$pb_quartiles) == FALSE, ], pb_quartiles), 
                    obs_ratio = length(atopally) / nrow(dat[is.na(dat$pb_quartiles) == FALSE, ]), 
                    obs_count = length(atopally),
                    democracies_mean = mean(jointpol7, na.rm = TRUE), 
                    rivals_mean = mean(rivalry_th, na.rm = TRUE), 
                    idealpointdiff_mean = mean(idealpointdiff_std, na.rm = TRUE), 
                    allies_mean = mean(atopally, na.rm = TRUE), 
                    cincdif_mean = mean(cincdif_std, na.rm = TRUE))
names(pb_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential")
pb_sum$`IGO type` <- 4 # "Peace-brokering IGOs"

igo_sum <- rbind(hligo_sum, ingr23_sum, mp3_sum, pb_sum)
names(igo_sum) <- c("Joint IGO memberships", "% of observations", "Observations", "Joint democracy", "Strategic rivalry", "UNGA ideal point difference", "Allies", "Power differential", "IGO type")
igo_sum$`IGO type` <- factor(igo_sum$`IGO type`, labels = c("IGOs with high leverage", "Structured IGOs", "Highly structured Security IGOs", "Peace-brokering IGOs"))
igo_sum <- igo_sum[is.na(igo_sum$`Joint IGO memberships`) == FALSE, ]

igo_sum$`% of observations` <- NULL
igo_sum_long <- reshape2::melt(igo_sum, id.vars = c("Joint IGO memberships", "IGO type", "Observations"))

# igo_sum$`Observations` <- NULL
# igo_sum_long <- reshape2::melt(igo_sum, id.vars = c("Joint IGO memberships", "IGO type", "% of observations"))

p <- ggplot(data = igo_sum_long, aes(x = `Joint IGO memberships`, y = value)) + facet_grid(variable ~ `IGO type`) + geom_line(aes(group = 1)) + geom_point(aes(size = Observations), shape = 21, fill = "white") + theme_jk() + xlab("Count of joint IGO memberships (quartile-based bins)") + ylab("Average (proportions or means standardized to 0-1 range)") + labs(size = "Claims\nin bin")

ggsave(p, file = "Output_Tables-and-Figures/claims_igo-profiles.pdf", width = 9.75, height = 10)

#############################################
# 4. Regression estimates for claims, with multiple claims per dyad-year reduced to most violent claim
#############################################

# Model 1

m1_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat_nomult)

m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
                       data = (m1_mf),
                       # family = binomial(link = "logit"),
                       # prior = student_t(df = 7, location = 0, scale = 2.5), 
                       # prior_intercept = normal(0, 10),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/claims_nomult_hiviol_m1_stanglm.RDS")

# Table
m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:9, 1, 11), ])
m1_tab$varname <- c("IGOs with high leverage", 
                    "Intangible salience of claim", 
                    "Tangible salience of claim", 
                    "Territorial claim", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/claims_nomult_hiviol_m1_table.tex", include.rownames = FALSE)

# First differences
m1_mm <- model.matrix(m1_eq, data = dat_nomult)

temp_sims <- as.matrix(m1_stanglm)
# temp_sims <- temp_sims[, grep("alpha|beta", colnames(temp_sims))]
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                          credint = c(0.025, 0.975), 
                          percentiles = c(0.1, 0.9),
                          full_sims = FALSE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]   [,2]
# [1,]     NA     NA
# [2,] 0.0000 4.0000
# [3,] 1.0000 2.3000
# [4,] 1.0000 5.0000
# [5,] 0.0000 1.0000
# [6,] 0.0000 1.0000
# [7,] 0.0000 1.0000
# [8,] 0.1496 2.6827
# [9,] 0.0000 1.0000

fd.dat <- data.frame(m1.fd)
fd.dat$variable_label <- c("IGOs with high leverage\n(from 0 to 4)", 
                           "Intangible salience of claim\n(from 1 to 2.3)", 
                           "Tangible salience of claim\n(from 1 to 5)", 
                           "Territorial claim\n(no to yes)", 
                           "Joint democracy\n(no to yes)", 
                           "Strategic rivalry\n(no to yes)",
                           "UNGA ideal point difference\n(0.1 to 2.7)",
                           "Allies \n(no to yes)")

# Alternative: plot posterior density

m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                          credint = c(0.025, 0.975), 
                          percentiles = c(0.1, 0.9),
                          full_sims = TRUE)

fd_long <- melt(m1.fd)
fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", 1,
                         ifelse(fd_long$Var2 == "salint", 2, 
                                ifelse(fd_long$Var2 == "saltan", 3,
                                       ifelse(fd_long$Var2 == "terriss", 4,
                                              ifelse(fd_long$Var2 == "jointpol7", 5,
                                                     ifelse(fd_long$Var2 == "rivalry_th", 6, 
                                                            ifelse(fd_long$Var2 == "absidealdiff", 7,
                                                                   ifelse(fd_long$Var2 == "atopally", 8, NA))))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", "IGOs with high leverage\n(from 0 to 4)",
                                 ifelse(fd_long$Var2 == "salint", "Intangible salience of claim\n(from 1 to 2.3)", 
                                        ifelse(fd_long$Var2 == "saltan", "Tangible salience of claim\n(from 1 to 5)",
                                               ifelse(fd_long$Var2 == "terriss", "Territorial claim\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)", 
                                                                    ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 2.7)",
                                                                           ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)", NA))))))))

# ROPE: use SE of the ratio of 1s

r_n <- length(m1_mf$useforce)
r_p <- sum(m1_mf$useforce) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  # scale_x_discrete(labels = rev(fd.dat$variable_label)) +
  # geom_rect(aes(ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf), col = "gray") +
  # geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of using force\nduring a claim as each predictor changes as indicated") + 
  xlab("") + 
  # scale_y_discrete(labels = rev(fd.dat$variable_label), expand = c(0, 0.6)) +
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Use of force less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Use of force more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.8, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-3, 3]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = ROPE, x = 7.25, xend = 7.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 8.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 8.25, xend = 8.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 3.5, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 3.45, xend = 3.15, size = 0.25) +
  annotate(geom = "text", y = -0.5, x = 7.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.275, x = 7.65, xend = 7.9, size = 0.25)

# % < 0

m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.1, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p, file = "Output_Tables-and-Figures/claims_nomult_hiviol_m1_fd.pdf", width = 6, height = 8)

# Compare LOO
m1_noigo <- stan_glm(formula = useforce ~ salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally,
                     data = (m1_mf),
                     family = binomial(link = "logit"),
                     prior = normal(0, 10), 
                     prior_intercept = normal(0, 10),
                     seed = 20170207,
                     cores = 4)

loo_noigo <- loo(m1_noigo, cores = 2)
loo_m1 <- loo(m1_stanglm, cores = 2)
compare(loo_noigo, loo_m1)
# elpd_diff se 
# 1.2       2.3 

# Influence of ind'l observations
p <- ggplot(data = data.frame(x = 1:151, y = loo_m1$pareto_k),
            aes(x = x, y = y)) + geom_point(shape = 1) + xlab("Observation") + ylab("Shape parameter k") + geom_hline(yintercept = 0.5, linetype = "dashed") + theme_jk()
ggsave(p, file = "Output_Tables-and-Figures/claims_nomult_hiviol_loo_k.pdf", width = 5, height = 4)

# Other model fit info with Stan

# Data 
roc_mf <-  model.frame(useforce ~ igo_lev3_count_use_bc + igo_ingr23_use_bc + igo_mp3_use_bc + igo_pb_use + igo_allothers_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + name + claim, data = dat_nomult)
N <- nrow(roc_mf)
useforce <- roc_mf$useforce
igo_lev3_count_use_bc <- roc_mf$igo_lev3_count_use_bc
igo_ingr23_use_bc <- roc_mf$igo_ingr23_use_bc
igo_mp3_use_bc <- roc_mf$igo_mp3_use_bc
igo_pb_use <- roc_mf$igo_pb_use
igo_allothers_bc <- roc_mf$igo_allothers_bc
salint <- roc_mf$salint
saltan <- roc_mf$saltan
terriss <- roc_mf$terriss
jointpol7 <- roc_mf$jointpol7
rivalry_th <- roc_mf$rivalry_th
absidealdiff <- roc_mf$absidealdiff
atopally <- roc_mf$atopally

m1_roc_data_list <- c("N", "useforce", "igo_lev3_count_use_bc", "igo_ingr23_use_bc", "igo_mp3_use_bc", "igo_pb_use", "igo_allothers_bc", "salint", "saltan", "terriss", "jointpol7", "rivalry_th", "absidealdiff", "atopally")

# Fit m1 directly in Stan
m1_roc_rstan <- stan(file = "claims_m1_roc.stan", 
                     data = m1_roc_data_list,
                     iter = 1000, chains = 4,
                     cores = 4)

print(m1_roc_rstan)

# Fit igo_ingr23_use_bc directly in Stan
m_ingr23_roc_rstan <- stan(file = "claims_m-ingr23_roc.stan", 
                           data = m1_roc_data_list,
                           iter = 1000, chains = 4,
                           cores = 4)

# Fit igo_mp3_use_bc directly in Stan
m_mp3_roc_rstan <- stan(file = "claims_m-mp3_roc.stan", 
                        data = m1_roc_data_list,
                        iter = 1000, chains = 4,
                        cores = 4)

# Fit igo_pb_use directly in Stan
m_pb_roc_rstan <- stan(file = "claims_m-pb_roc.stan", 
                       data = m1_roc_data_list,
                       iter = 1000, chains = 4,
                       cores = 4)

# Fit all other IGOs directly in Stan
m_allothers_roc_rstan <- stan(file = "claims_m-allothers_roc.stan", 
                              data = m1_roc_data_list,
                              iter = 1000, chains = 4,
                              cores = 4)

# Fit empty model (no IGOs) directly in Stan
m_noigos_roc_rstan <- stan(file = "claims_m-noigos_roc.stan", 
                           data = m1_roc_data_list,
                           iter = 1000, chains = 4)

# ROC and precision-recall curves

m1_rp <- MCMC_roc_prc(stan_object = m1_roc_rstan, 
                      model_frame = m1_mf, 
                      model_name = "IGOs with high leverage")

m_ingr23_rp <- MCMC_roc_prc(stan_object = m_ingr23_roc_rstan, 
                            model_frame = m1_mf, 
                            model_name = "Structured IGOs")

m_mp3_rp <- MCMC_roc_prc(stan_object = m_mp3_roc_rstan, 
                         model_frame = m1_mf, 
                         model_name = "Highly structured Security IGOs")

m_pb_rp <- MCMC_roc_prc(stan_object = m_pb_roc_rstan, 
                        model_frame = m1_mf, 
                        model_name = "Peace-brokering IGOs")

m_allothers_rp <- MCMC_roc_prc(stan_object = m_allothers_roc_rstan, 
                               model_frame = m1_mf, 
                               model_name = "All other IGOs")

m_noigos_rp <- MCMC_roc_prc(stan_object = m_noigos_roc_rstan, 
                            model_frame = m1_mf, 
                            model_name = "No IGOs")

PRC_dat <- rbind(m1_rp$prc_dat, m_ingr23_rp$prc_dat, m_mp3_rp$prc_dat, m_pb_rp$prc_dat, m_allothers_rp$prc_dat, m_noigos_rp$prc_dat)
PRC_labeldat <- rbind(m1_rp$area_under_prc, m_ingr23_rp$area_under_prc, m_mp3_rp$area_under_prc, m_pb_rp$area_under_prc, m_allothers_rp$area_under_prc, m_noigos_rp$area_under_prc)

PRC.p <- ggplot(data = PRC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = PRC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 0, y = 0.2, hjust = 0) + labs(title = "Precision-recall curves", hjust = 0.5) + xlab("Recall") + ylab("Precision") + theme_jk()

ROC_dat <- rbind(m1_rp$roc_dat, m_ingr23_rp$roc_dat, m_mp3_rp$roc_dat, m_pb_rp$roc_dat, m_allothers_rp$roc_dat, m_noigos_rp$roc_dat)
ROC_labeldat <- rbind(m1_rp$area_under_roc, m_ingr23_rp$area_under_roc, m_mp3_rp$area_under_roc, m_pb_rp$area_under_roc, m_allothers_rp$area_under_roc, m_noigos_rp$area_under_roc)

ROC.p <- ggplot(data = ROC_dat, aes(x = x, y = y)) + geom_line() + ylim(c(0, 1)) + facet_wrap(~ Model, ncol = 1) + geom_text(data = ROC_labeldat, aes(label = paste("Area under\ncurve:",round(auc, digits = 2))), x = 1, y = 0.2, hjust = 1) + labs(title = "ROC curves", hjust = 0.5) + xlab("1 - Specificity") + ylab("Sensitivity") + theme_jk()

pdf(file = "Output_Tables-and-Figures/claims_nomult_hiviol_pr-roc.pdf", width = 6, height = 12)
grid.arrange(PRC.p, ROC.p, ncol = 2)
dev.off()

# Posterior predictive check

stan_object <- m1_roc_rstan
pred_sim <- as.data.frame(stan_object, pars = c("pred"))
pred_sim <- reshape2::melt(pred_sim)
# pred_sim$obs <- as.numeric(gsub("[^0-9]","",as.character(pred_sim$Parameter)))
pred_sim$obs <- as.numeric(gsub("[^0-9]","",as.character(pred_sim$variable)))

obs_igo <- data.frame(y_obs = m1_mf[, "useforce"], 
                      x_obs = m1_mf[, "igo_lev3_count_use_bc"],
                      obs = 1:(nrow(m1_mf)))

median_obs <- summarize(group_by(obs_igo, x_obs),
                        median_y = mean(y_obs))

pp_df <- merge(x = pred_sim[, c("value", "obs")], y = obs_igo, by = "obs")

ppc.p <- ggplot(data = pp_df, aes(x = x_obs, y = value, group = x_obs)) + geom_boxplot(outlier.size = NA, size = 0.25) + scale_x_continuous(breaks = unique(obs_igo$x_obs)) + geom_point(data = median_obs, aes(x = x_obs, y = median_y), color = "red") + theme_jk() + xlab("Joint memberships in IGOs with high leverage") + ylab("Estimated probability of using force") + annotate("point", x = 0.2, y = 0.9, color = "red") + annotate("text", x = 0.4, y = 0.9, hjust = 0, label = "Percent of observed claims\nexperiencing the use of force", size = 3) + annotate("boxplot", x = 5, lower = 0.85, upper = 0.95, middle = 0.9, ymin = 0.8, ymax = 1, size = 0.25) + annotate("text", x = 5.5, y = 0.9, hjust = 0, label = "Posterior predictive\ndistribution", size = 3)

ggsave(ppc.p, file = "Output_Tables-and-Figures/claims_nomult_hiviol_ppc.pdf", width = 6, height = 4)

## BMA

bma_dat_nomult <- model.frame(useforce ~ igo_lev3_count_use_bc + igo_ingr23_use_bc + igo_mp3_use_bc + igo_pb_use + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally + cincdif, data = dat_nomult)
names(bma_dat_nomult) <- c("Use of force",
                           "IGOs with high leverage",
                           "Structured IGOs",
                           "Highly structured Security IGOs", 
                           "Peace-brokering IGOs",
                           "Intangible salience of claim", 
                           "Tangible salience of claim", 
                           "Territorial claim", 
                           "Joint democracy", 
                           "Strategic rivalry",
                           "UNGA ideal point difference",
                           "Allies",
                           "Power differential")

claims_nomult_bma <- bic.glm(y = bma_dat[, 1], x = bma_dat[, c(2:13)], 
                             glm.family = binomial(link = "logit"))

pdf("Output_Tables-and-Figures/claims_nomult_hiviol_bma.pdf", width = 12, height = 8)
plotJK.bic.glm(claims_nomult_bma, mfrow = c(3, 4))
dev.off()

# End of file.